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Abstract 


For  3D  surface  reconstruction  problems  with  noisy  and  incomplete  range  data  measured 
from  complex  scenes  with  arbitrary  topologies,  a  low-level  representation,  such  as  level 
set  surfaces,  is  used.  Such  surface  reconstruction  is  typically  accomplished  by  minimiz¬ 
ing  a  weighted  sum  of  data-model  discrepancy  and  model  smoothness  terms.  This  paper 
introduces  a  new  nonlinear  model  smoothness  term  for  surface  reconstruction  based  on 
variations  of  the  surface  normals.  A  direct  solution  requires  solving  a  fourth-order  partial 
differential  equation  (PDE),  which  is  very  difficult  with  conventional  numerical  techniques. 
Our  solution  is  based  on  processing  the  normals  separately  from  the  surface,  which  allows 
us  to  separate  the  problem  into  two  second-order  PDEs.  The  proposed  method  can  smooth 
complex,  noisy  surfaces,  while  preserving  sharp,  geometric  features,  and  it  is  a  natural  gen¬ 
eralization  of  edge -preserving  methods  in  image  processing,  such  as  anisotropic  diffusion. 


Chapter  1 
Introduction 


The  precision  of  range  measurement  systems,  such  as  time-of-flight  laser  range  finders, 
has  been  increasing  while  their  price  drops.  If  combined  with  a  well-founded  method  for 
surface  reconstruction,  these  improvements  could  make  capturing  3D  shape  as  ubiquitous 
as  photography.  However,  significant  challenges  to  surface  reconstruction,  such  as  mea¬ 
surement  noise  and  variations  in  measurement  density,  remain.  This  paper  addresses  the 
problem  of  preserving  geometric  features,  i.e.  edges,  corners  and  junctions  on  surfaces,  in 
full  3D  reconstructions  of  complex  scenes  from  multiple,  noisy  range  images.  Figure  1.1(a) 
illustrates  a  typical  range  image.  Measurement  noise  and  occlusions  (between  the  file  cab¬ 
inet  and  the  chair)  can  be  observed  in  this  data.  Figure  1.1(b)  illustrates  the  reconstruction 
with  the  proposed  method  from  a  similar  view  point.  The  result  is  an  improvement  over 
the  state-of-the-art  full  3D  reconstructions  for  complex  scenes:  creases  and  comers  at  the 
interchapters  of  the  various  planes  in  the  scene  have  been  preserved  while  measurement 
noise  has  been  effectively  eliminated. 


(a)  (b) 

Figure  1.1:  (a)  A  range  image  surface,  and  (b)  feature  preserving  surface  reconstruction. 
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Full  3D  reconstruction  recovers  a  view-independent  surface  model  from  multiple  registered 
range  images.  This  problem  is  distinct  from  depth  reconstruction ,  also  known  as  2|D 
reconstruction,  which  recovers  structure  from  only  one  point  of  view  or  a  stereo  pair  of 
images  [1,  2,  3].  Depth  reconstruction  does  not  produce  a  model  that  makes  sense  when 
viewed  from  different  viewpoints  or  when  determining  inherently  3D  properties,  such  as 
volume.  For  instance,  occluded  portions  of  the  scene  in  Figure  1.1(a)  are  present  in  the 
reconstructed  model  because  multiple  range  images  are  used  in  full  3D  reconstruction. 
This  result  would  not  have  been  possible  with  depth  reconstruction  methods.  The  full  3D 
problem  is  not  a  mere  extension  of  depth  reconstruction  because  it  lacks  the  following 
properties  of  the  latter:  the  depth  map  has  a  well-defined  topology  (a  function  of  two 
variables)  and  there  is  a  one-to-one  correspondence  between  the  measurements  and  the 
positions  on  the  model. 

Recovering  a  full  3D  model  from  a  set  of  noisy  2D  range  images  is  an  ill-posed  inverse 
problem.  Hence,  the  estimator  can  not  depend  solely  on  the  input  data,  and  requires  regu¬ 
larization.  Regularization  reduces  the  effects  of  measurement  noise  and  fills  surfaces  in  a 
plausible  manner  where  there  is  no  data  from  any  of  the  range  images  by  placing  additional 
constraints  on  the  reconstructed  surfaces.  This  problem  has  been  approached  in  the  com¬ 
puter  vision  literature  mainly  as  a  problem  of  finding  sets  of  geometric  primitives  that  best 
represent  the  objects  being  measured  [4,  5,  6,  7].  Primitives  typically  have  only  a  few  shape 
parameters,  i.e.,  height  and  radius  for  a  cylinder;  therefore,  impose  their  own  structure  on 
to  the  data.  In  this  case  regularization  is  inherent  to  the  surface  model.  Such  approaches  are 
suitable  for  higher- level  tasks  of  object  recognition  and  decomposition  into  parts;  however, 
they  are  limited  to  modeling  relatively  simple  objects. 

An  alternative  is  to  use  level  set  surfaces  [8],  which  are  a  non-parametric  shape  represen¬ 
tation.  Level  set  surfaces  can  be  used  to  recover  scenes  with  arbitrarily  complex  geometry 
and  topology.  The  reconstructed  level  set  models  are  not  limited  to  prescribed  topologies, 
and  can  adapt  to  the  topology  of  the  measured  scenes  automatically.  However,  level  set  sur¬ 
faces  do  not  have  a  rigid  shape  structure,  and  therefore,  regularization  must  be  performed 
explicitly. 

We  formulate  surface  estimation  (reconstruction)  in  a  variational  energy  optimization  frame¬ 
work.  Variational  methods  typically  minimize  an  energy  function,  which  is  a  weighted  sum 
of  an  input  data-model  discrepancy  term  and  a  model  smoothness  term,  with  respect  to  the 
model.  Then  the  surface  estimator  is  defined  as 

S(D,  a)  =  arginf  [F(S,  D)  +  aP(S )]  (1.1) 

where  a  determines  the  relative  weights  of  the  terms.  The  input  data-model  discrepancy 
term,  F(S,  D ),  forces  the  surface  estimator  to  be  “close”  to  the  measured  data.  The  model 
smoothness  term,  P(S),  provides  regularization. 


This  paper  studies  the  model  smoothness  term.  Surface  area  is  a  simple  choice  and  has  been 
used  extensively  in  previous  work  [9,  10].  This  is  based  on  the  assumption  that  among 
surfaces  with  similar  data-model  discrepancy  measures,  those  that  have  smaller  area  are 
simpler  than  surfaces  of  larger  area,  and  therefore  have  a  higher  chance  of  occurrence  in 
reality.  Despite  its  simplicity,  surface  area  as  a  measure  of  smoothness  has  significant 
drawbacks,  such  as  pinching  of  thin  structures.  We  argue  that  measuring  the  variation 
of  the  surface  normal  vectors  offers  a  better  and  more  flexible  alternative.  Specifically, 
this  family  of  smoothness  functions  offers  an  elegant  and  mathematically  correct  general¬ 
ization  of  Perona  &  Malik’s  (P&M)  anisotropic  diffusion  [12]  to  surface  reconstruction. 
This  generalization  allows  us  to  preserve  important  shape  features  such  as  edges,  corners 
and  junctions  in  reconstructed  scenes  while  effectively  eliminating  measurement  noise  and 
other  artifacts. 

The  remainder  of  this  paper  is  organized  as  follows.  Chapter  2  discusses  the  related  work 
in  the  literature.  Chapter  3  discusses  level  set  surface  reconstruction  and  proposes  a  gener¬ 
alization  of  P&M’s  edge  detection  method  as  a  level  set  surface  regularization  term.  Chap¬ 
ter  4  solves  the  proposed  regularization  term  as  a  level  set  motion.  Chapter  5  demonstrates 
the  quantitative  advantages  of  the  proposed  method  and  provides  examples  of  reconstruc¬ 
tion  of  real,  complex  scenes  from  noisy  range  data.  Chapter  6  summarizes  the  contributions 
of  this  paper  and  discusses  possibilities  for  future  research  directions. 


Chapter  2 
Related  Work 


As  mentioned  earlier,  computer  vision  researchers  approach  surface  reconstruction  either 
as  a  depth  reconstruction  problem  [1,  2,  3]  or  a  view  independent  problem.  Earlier  literature 
on  the  view  independent  problem  focuses  mainly  on  high-level  approaches  that  fit  various 
geometric  primitives  to  the  data  [4,  5,  6,  7].  Both  problems  are  different  from  a  view 
independent  reconstruction  of  complex  scenes  with  low-level  shape  representations.  In 
response  to  the  advances  in  3D  range  sensing  device,  researchers  in  a  variety  of  fields 
have  started  to  study  this  problem.  In  computer  graphics,  the  accuracy  of  the  data  exceeds 
the  requirements  of  the  application,  and  therefore  the  problem  is  treated  as  a  problem  of 
assembling  pieces  of  noiseless  information.  For  instance,  Turk  and  Levoy  [13]  propose 
a  “zippering”  algorithm  for  combining  triangle  meshes  of  range  maps  of  an  object  from 
different  points  of  view.  Curless  and  Levoy  [14]  take  into  account  measurement  noise  by 
averaging  range  information  in  a  volumetric  representation.  However,  their  method  is  not 
based  on  statistics  the  scanner  and  model  geometry. 

Several  authors  [9,  10]  demonstrate  the  advantages  of  using  level  set  methods  for  recon¬ 
structing  complex  shapes.  They  use  mean  curvature  flow,  a  second-order  level  set  partial 
differential  equation  (PDE),  to  obtain  a  smooth  solution.  This  flow  is  the  gradient  de¬ 
scent  for  the  first  variation  of  surface  area  [15].  Hence,  for  regularization  purposes,  both 
approaches  are  formulated  with  the  surface  area  model  smoothness  term  in  the  estimator 
(1.1).  Mean  curvature  flow  suffers  from  several  problems  including  volume  shrinkage, 
pinching  of  thin  structures,  and  elimination  of  sharp  features  [11].  These  problems  can  be 
observed  in  the  results  presented  in  Chapter  5. 

To  alleviate  the  problems  associated  with  mean  curvature  flow,  several  authors  have  pro¬ 
posed  smoothing  level  set  surfaces  by  modified  second-order  flows  that  use  weighted  com- 
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binations  of  principle  curvatures.  For  instance,  Lorigo  et  al.  [16]  propose  a  smoothing 
flow  that  uses  the  minimum  curvature  for  tubular  structures.  Clarenz  et  al.  [17]  propose  an 
anisotropic  surface  mesh  diffusion  as  a  modified  second  order  flow,  but  this  modified  PDE 
lacks  a  variational  basis,  and  therefore  is  not  useful  for  surface  reconstruction.  Steven¬ 
son  and  Delp  propose  minimizing  total  curvature,  which  is  the  surface  integral  of  the  sum 
of  the  squared  principal  curvatures,  for  regularizing  depth  reconstruction  and  parametric 
Monge  patches  [18].  This  method  does  not  apply  to  level  set  surfaces.  Furthermore,  the 
gradient  descent  flow  for  total  curvature  is  a  fourth-order  PDE  that  is  computationally  ex¬ 
pensive  and  unstable  to  compute.  Hence,  Stevenson  and  Delp  use  a  non-geometric  thin 
plate  approximation. 

In  previous  work  [11],  we  propose  a  two-step  approach  to  surface  smoothing:  (1)  operate 
on  the  normal  map  of  a  surface,  and  (2)  refit  a  surface  to  the  processed  normals.  In  this 
paper,  we  show  that  an  quadratic  measure  on  the  variations  in  surface  normals  variations 
is  equivalent  to  total  curvature.  Unlike  using  calculus  of  variations  on  the  total  curvature 
of  the  surface  directly,  our  formulation  in  terms  of  the  normals  yields  a  computationally 
tractable  and  stable  gradient  descent  flow.  Furthermore,  we  also  propose  a  robust  measure 
of  surface  normal  variations  that  allows  a  novel  generalization  of  P&M  feature  preserving 
anisotropic  image  diffusion  [12]  to  surface  reconstruction. 

Our  energy  optimization  approach  can  also  be  stated  in  terms  of  Bayesian  maximum  a 
posteriori  (MAP)  estimation.  According  to  Bayes  rule,  MAP  estimators  maximize  the 
logarithm  of  the  product  of  two  distinct  probabilities:  the  likelihood  of  the  measurement 
data  conditioned  on  the  surface  model  and  the  prior  probability  distribution  of  the  model. 
In  the  estimator  given  by  (1.1),  the  input  data-model  discrepancy  term  corresponds  to  the 
logarithm  of  the  conditional  likelihood  and  the  model  smoothness  term  corresponds  to  the 
logarithm  of  the  prior  [10].  This  brings  up  a  connection  with  other  approaches  that  use 
shape  priors  for  active  contour  and  level  set  models  [19,  20,  21,  22].  However,  in  all  of 
these  works,  the  “prior”  is  a  global  description  of  expected  shape(s)  learned  from  a  training 
set.  In  the  segmentation  and/or  registration  stage,  the  shape  model  is  forced  to  be  a  rigid 
transformation  of  the  learned  prior  shape  with  some  tolerance  for  local  variations.  This 
approach  is  useful  for  reconstructing/segmenting  specific  classes  of  shapes,  such  as  cortical 
surfaces  from  head  MRI  data;  however,  learned  priors  can  not  be  used  for  regularization 
of  reconstructed  models  in  a  general  setting.  The  commonality  between  general  shapes 
is  not  on  a  global  level,  but  on  a  lower  level,  such  as  probability  distributions  for  surface 
normal  variations,  which  we  use  in  this  paper.  The  quadratic  and  robust  penalty  term  on 
surface  normal  variations  yield  generic  isotropic  and  anisotropic  low  level  shape  priors, 
respectively.  These  priors  are  not  learned  from  a  training  set. 


Chapter  3 

Variational  Implicit  Surface 
Reconstruction 


A  deformable  surface,  S(t),  can  be  represented  as  the  zero  level  set  of  a  higher  dimensional 
embedding  function  (j)  :  R3  x  R  ,  R,  S(t)  =  |x(f)  €  R3  ]  <j)  (x(t).  f)  =  0},  where  t  is 
the  evolution  parameter  (time).  Surfaces  defined  in  this  way  divide  a  volume  into  two  parts: 
inside  (</)  >  0)  and  outside  (<p  <  0).  The  family  of  PDEs  that  describe  motions  of  S  via 
d(j)/dt,  and  the  upwind  scheme  for  solving  them  on  a  discrete  grid  is  the  methods  of  level 
sets  [8]. 

The  surface  reconstruction  energy  of  (1.1)  can  be  expressed  as  a  function  of  the  data  D, 
and  the  level  set  model  (f).  Whitaker  [10]  formulates  the  input  data-model  discrepancy  part 
of  this  energy  as  a  volumetric  integral 

F(<t>,  D)  =  [  G(x,£>)ff(^(x))dx,  (3.1) 

Jfi 

where  f)  is  the  volumetric  domain  swept  out  by  the  range  images,  G  is  an  accurate  line- 
of-sight  error  function,  and  H  is  the  Heaviside  function  [23].  Minimizing  (3.1)  by  itself 
would  correspond  to  a  maximum  likelihood  estimator ,  but  a  maximum  apriori  estimator  is 
implemented  by  using  surface  area  as  a  model  smoothness  term: 

P(cf))=  [  ||V0(x)||  dx.  (3.2) 

J  n 

Then,  the  gradient  descent  for  (j)  is 


dcf)/dt  =  ||V</>(f)J|  (G(x,  D)  +  aif(x)) , 
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(3.3) 


(a)  (b) 

Figure  3.1:  (a)  A  quadratic  penalty  term  with  a  hard  cutoff,  (b)  an  exponential  penalty  term 
corresponding  to  P&M’s  anisotropic  diffusion. 

where  H  =  V  •  (V</>/  j  Voj  |)  is  mean  curvature  [15].  Using  discrete  time  steps,  the  model 
is  evolved  as  cf)(t  +  At)  =  <p(t)  +  At  d(f)/dt.  The  steady-state  solution  of  this  evolution  is 
the  surface  estimator:  S  =  {x  e  R3  |  lim^oo  <p  (x,  t)  =  0}. 

In  this  paper,  we  use  the  same  data-model  discrepancy  term  and  the  initialization  method 
for  cj){t  =  0).  However,  we  propose  a  better  model  smoothness  term  based  on  measures  on 
the  variation  of  the  surface  normal  vectors,  N, 

P(4 0  =  J  E  (l|V*N(x)||/r)  ||V0(x)||  Ac,  (3.4) 

where  V^N  is  the  matrix  whose  rows  are  the  gradient  vectors  of  the  respective  components 
of  N  intrinsic  to  the  isosurfaces  of  <j>.  The  Frobenius  matrix  norm,  i.e.  the  square  root  of 
the  sum  of  squares  of  the  matrix  elements,  is  denoted  by  1 1  - 1 1  .  Notice  that  if  E(y)  =  1, 
(3.4)  reduces  to  surface  area  (3.2).  On  the  other  hand,  if  we  choose  a  quadratic  measure, 
E(y)  =  y2,  we  obtain  the  sum  of  squared  principal  curvatures  (total  curvature)  [11].  The 
resulting  gradient  descent  PDE  is  analogous  to  running  the  heat  equation  PDE  on  the  sur¬ 
face  normals.  Figure  4.1  (a)  illustrates  a  noisy  field  of  unit  vectors  which  are  the  gradient 
directions  for  the  noisy  distance  transform  of  the  polygon  shown.  Figure  4.1  (b)  demon¬ 
strates  the  result  of  smoothing  with  the  PDE  derived  above.  Similar  to  the  heat  equation, 
this  flow  eliminates  both  noise  and  sharp  discontinuities  in  the  data.  Our  goal  is  to  preserve 
these  discontinuities  while  penalizing  noisy  variations  of  the  normals  elsewhere.  This  is 
similar  to  segmenting  the  normal  map. 

Mumford  and  Shah  formulate  the  problem  of  image  segmentation  in  a  variational  frame¬ 
work  [24].  The  Mumford-Shah  energy  is  the  sum  of  three  terms:  (i)  the  data-model  dis¬ 
crepancy,  (ii)  quadratic  penalty  on  model  image  smoothness  over  the  image  domain  except 
on  a  set  of  discontinuities  modeled  by  a  binary  image,  and  (iii)  the  length  of  the  discontinu¬ 
ities  in  that  binary  image.  The  sum  of  the  latter  two  terms  correspond  to  using  the  penalty 
function,  E(y),  shown  in  Figure  3.1(a).  The  existence  of  a  binary  model  poses  serious  dif- 


Acuities  in  the  estimation  process.  Nordstrom  [25]  and  Black  et  al.  [26]  show  that  P&M 
anisotropic  diffusion  approach  to  edge  detection  [12]  is  equivalent  to  using  a  corresponding 
robust  penalty  term  in  the  Mumford-Shah  segmentation  framework.  This  penalty  term, 

E(y)  =  l-e~^,  (3.5) 

shown  in  Figure  3.1(b),  avoids  using  a  binary  discontinuity  model.  The  parameter  /i  con¬ 
trols  the  degree  of  edge  preservation.  The  above  result  can  readily  be  generalized  to  normal 
vectors  by  substituting  y  =  ||V^N(x)||^r.  Figure  4.1(c)  illustrates  the  result  of  smooth¬ 
ing  the  noisy  image  of  normals  with  this  penalty  term.  The  discontinuities  in  the  normal 
directions  between  the  quadrants  are  preserved  while  the  noise  is  smoothed.  Minimization 
of  energies  of  the  form  (3.4)  require  solving  fourth-order  level  set  PDEs,  a  computation¬ 
ally  unstable  and  expensive  task.  The  next  chapter  introduces  a  method  for  breaking  the 
solution  into  two  second-order  PDEs  that  can  be  efficiently  solved. 


Chapter  4 

Level  set  motion  via  normal  map 
diffusion 


In  equation  (3.4),  Y?)N  is  the  matrix  whose  rows  are  the  gradient  vectors  of  the  respective 
components  of  N  intrinsic  to  the  isosurfaces  of  (j).  By  intrinsic,  we  mean  that  when  using 
implicit  representations  one  must  account  for  the  fact  that  derivatives  of  functions  defined 
on  the  surface  are  computed  by  projecting  their  3D  derivatives  onto  the  surface  tangent 
plane.  The  3  x3  projection  matrix  for  the  implicit  surface  normal  is  P  =  \V<p\  2, 

where  ®  is  the  tensor  product.  The  projection  matrix  onto  the  surface  tangent  plane  is  I  — P, 
where  I  is  the  identity  matrix.  Then  the  intrinsic  gradient  of  the  normals  can  be  defined 
using  this  projection  operator  and  the  regular  Euclidean  gradient  V^N  =  VN  (I  —  P). 

Given  an  initialization  for  (j)  with  the  methods  described  in  [10],  we  compute  the  surface 
normals  N  =  V</>/  ||V</>||.  Then,  to  avoid  solving  fourth-order  level  set  PDEs  directly, 


(a) 


(b) 


iv.  ^  ^  ^  tv.  St 


Figure  4.1:  (a)Normals  computed  from  a  noisy  distance  transform  of  a  polygon,  normals 
smoothed  with  (b)  the  quadratic  penalty  term,  and  (c)  the  robust  penalty  term. 
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we  decouple  N  from  <j>.  In  other  words,  we  fix  <p  (the  surface  shape)  as  we  process  the 
normals  to  minimize  the  energy  given  by  (3.4).  Solutions  to  constrained  minimization  for 
unit  vectors  on  an  implicit  surface  are  discussed  in  [27];  however,  the  goal  in  that  work  is  to 
not  to  smooth  surfaces,  but  to  diffuse  general  vector  functions  on  surfaces.  We  implement 
the  constrained  minimization  with  the  following  second-order  PDE: 


dN 

dt 


(I 


N  <g>  N)  V 


E' 


(4.1) 


where  E'(y)  is  the  derivative  of  E(y).  For  the  penalty  term  given  in  (3.5),  g  =  e  i*7  up 
to  a  constant  multiplicative  factor.  Figure  4.1  (a)  illustrates  a  noisy  field  of  unit  vectors. 
Figure  4.1  (b)  and  (c)  demonstrate  results  of  smoothing  with  the  choices  of  E(y)  discussed 
in  Chapter  3. 


The  next  step  is  to  relate  the  deformation  of  the  level  sets  of  (j)  to  the  evolution  of  N.  Once 
more,  using  a  variational  approach,  we  can  manipulate  (t>  so  that  it  fits  the  new  normal  field 
by  minimizing  a  penalty  function, 


v{4>)=  [  L/V0-  W-  v<0-  vn|  dS , 
Jn  L  J 


(4.2) 


that  quantifies  the  discrepancy  between  the  gradient  vectors  of  <j)  and  the  target  normal  map 
B allester  et  al.  [28]  use  the  same  function  for  filling  in  missing  regions  in  images  by  joint 
interpolation  of  the  image  intensity  and  its  gradient.  The  first  variation  of  this  function  with 
respect  to  <p  is 


d<fi 


V(j) 

LIIWI 


-  N 


-  [H*  -  tfN' 


(4.3) 


where  H®  is  the  mean  curvature  of  the  level  set  surface  and  HN  is  half  the  divergence  of 
the  normal  map.  Figure  4.1(b)  and  (c)  illustrate  the  curves  refitted  to  the  smoothed  normal 
fields  with  this  approach.  Finally,  the  gradient  descent  for  the  surface  reconstruction  with 
the  model  smoothness  energy  (3.4)  is 


W/dt  =  ||V<Mt)||  (G(x,D)+a(H*(x)  -  tfN(x))) ,  (4.4) 


which  is  similar  to  (3.3),  but  has  a  different  smoothing  term. 


The  flow  chart  for  the  algorithm  is  shown  in  Figure  4.2.  We  have  derived  a  gradient  descent 
for  the  normal  map  that  minimize  the  energy  functions  of  the  form  (3.4).  The  normals 
processing  stage  of  the  algorithm  computes  the  gradient  descent  for  the  normals  defined  in 
(4.1)  for  a  fixed  number  of  iterations  (25  for  the  experiments  in  this  paper).  Hence,  we  avoid 
evolving  evolving  the  normals  too  far  away  from  their  initialization  from  (I).  The  surface 
fitting  to  the  the  combined  normal  map  and  data  terms  is  given  as  a  gradient  descent  in  (4.4). 


Figure  4.2:  Surface  reconstruction  flow  chart. 


This  stage  of  the  algorithm  is  run  until  the  discrepancy  measure  (4.2)  between  the  new 
normals  and  0  ceases  to  decrease,  which  signals  the  need  for  another  round  of  processing 
the  normal  vectors.  The  overall  algorithm  repeats  these  two  steps  to  minimize  the  surface 
reconstruction  energy  in  terms  of  (f>  until  the  RMS  value  for  d<p/dt  becomes  small  (less 
than  1CT6),  which  signals  convergence.  This  algorithm  consists  of  solving  two  second- 
order  PDEs  in  series  instead  of  a  direct  fourth-order  PDE,  which  makes  it  computationally 
tractable.  We  show  the  relationship  of  this  algorithm  to  solving  the  direct  fourth-order  PDE 

in  [11]. 


Chapter  5 
Experiments 


In  this  chapter,  we  compare  reconstructions  with  proposed  the  model  smoothness  energies 
against  reconstructions  with  the  standard  surface  area  energy.  For  the  proposed  family  of 

_y2 

energies  (3.4),  we  will  call  the  choice  of  E(y)  =  y2  and  E(y)  =  1  —  e  U,  the  isotropic  and 
anisotropic  reconstructions,  respectively.  Note  that  /j  is  fixed  at  0.2  for  all  the  experiments. 
Unlike,  in  P&M  image  diffusion,  this  parameter  does  not  need  to  be  changed  for  different 
surface  reconstructions.  In  the  context  of  P&M  image  diffusion,  the  units  of  jj,  are  in  gray 
levels;  consequently,  the  optimal  choice  of  /i  is  image  dependent.  In  surface  reconstruction, 
the  units  are  in  curvature  which  is  data  independent.  This  makes  it  possible  to  choose  a  /j 
value  that  gives  consistent  results  over  a  broad  range  of  surfaces. 

We  first  experiment  with  geometric  shapes  for  which  we  can  construct  analytical  distance 
transforms.  We  use  the  following  experiment  setup: 


1.  Build  range  images  by  simulating  the  laser  range  finder  located  at  several  positions, 

2.  Add  independent  Gaussian  noise  to  the  range  images, 

3.  Reconstruct  a  model,  and 

4.  Compute  the  root  mean  square  geometric  distance,  between  the  resulting  model  and 
the  analytical  shape. 


The  first  shape  we  examine  is  a  sphere  with  radius  1  unit.  All  other  distances  are  relative 
to  this  measurement  unit.  For  this  experiment  we  simulate  six  range  finders  located  at  a 
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distance  of  3.5  units  from  the  center  of  the  sphere  along  the  six  cardinal  directions.  In¬ 
dependent  Gaussian  noise  with  a  standard  deviation  of  0.1  units,  is  added  to  each  range 
image.  Figure  5.1(a)  plots  the  RMS  error,  E,  against  the  logarithm  of  the  regularization 
weight,  logo;,  for  the  different  reconstructions.  The  units  on  the  y-axis  are  the  same  as  the 
units  used  to  described  the  size  of  the  shape.  It  can  be  observed  from  Figure  5.1  that  the 
limiting  value  for  E  as  a  — >  0  is  approximately  0.0125.  This  limit  is  the  error  obtained  if 
surface  reconstruction  is  performed  without  regularization.  This  error  level  is  smaller  than 
the  noise  added  to  the  range  images  because  of  the  averaging  effect  of  using  multiple  range 
images.  The  anisotropic  and  the  isotropic  curvature  priors  at  their  optimal  weight  provide 
a  75%  reduction  on  this  error.  On  the  other  hand,  surface  area  provides  slightly  better  than 
a  50%  reduction  at  its  optimal  weight.  The  shapes  of  the  error  plots  is  more  important  than 
the  results  at  optimal  choices  of  weight.  The  surface  area  prior  performs  poorly  as  a  is 
increased  beyond  1;  this  is  due  to  the  fact  that  the  surface  area  prior  causes  shrinkage  in  the 
surface  models.  In  practice,  this  will  mean  difficulties  for  the  user  in  choosing  a  weight  for 
surface  area  regularization  that  works  different  reconstruction  scenarios.  In  contrast,  the 
proposed  reconstructions  have  relatively  flat  error  plots.  Isotropic  reconstruction  is  as  good 
as  the  anisotropic  reconstruction  because  the  sphere  does  not  contain  creases. 

To  examine  the  differences  between  isotropic  and  anisotropic  reconstruction  further,  we 
experiment  with  a  cube  shape.  In  this  experiment,  we  use  8  range  finder  locations  (one 
in  each  octant).  Figure  5.2  (a)  and  (b)  show  the  original  cube  with  sides  1  unit  long,  and 
the  surface  initialization  from  the  noisy  range  images,  respectively.  Independent  Gaussian 
noise  with  standard  deviation  0.1  was  added  to  the  simulated  data  to  create  the  noisy  range 
images.  The  results  (see  Figure  5.2)  with  isotropic  reconstruction  have  rounded  corners 
compared  to  the  successfully  denoised,  approximately  piecewise  planar  results  obtained 
with  the  anisotropic  reconstruction. 

The  next  example  involves  12  real  range  scans  of  a  room  which  were  registered  using  the 
methods  described  in  [10].  A  close-up  view  of  a  portion  of  one  of  the  range  images  and  the 
result  of  anisotropic  reconstruction  are  shown  in  Figure  1.1.  The  anisotropic  reconstruction 
of  the  entire  scene  is  shown  in  Figure  5.3.  We  now  examine  reconstructions  of  one  of  the 
chairs  in  this  scene.  Figure  5.4  (a),  (c)  and  (e)  illustrate  the  results  obtained  by  qualitatively 
choosing  good  values  for  a.  Figure  5.4(b),  (d)  and  (f)  illustrate  the  results  if  a  is  chosen 
to  be  10  times  this  value.  These  results  show  that  the  anisotropic  reconstruction  produces 
the  best  results  and  is  least  sensitive  to  the  choice  of  a.  Another  well  known  problem  with 
surface  area  reconstruction  can  easily  be  observed  in  Figure  5.4(b);  the  beam  connecting 
the  base  to  the  seat  is  being  pinched-off.  This  experiment  illustrates  the  importance  of  the 
anisotropic  reconstruction  for  scenes  with  high  curvature  features  and  sharp  creases. 

Figure  5.5  illustrates  anisotropic  reconstructions  of  a  vehicle.  This  example  further  illus¬ 
trates  the  success  of  the  anisotropic  reconstruction  in  denoising  data  with  sharp  features. 


-  -x-  -  Surface  area 
— • —  Isotropic  curvature 


-  Anisotropic  curvature 


(a)" 


Figure  5.1:  Rms  distance  between  the  reconstructed  and  the  analytical  surface  for  (a)  the 
sphere,  and  (b)  the  cube. 


Figure  5.2:  (a)  Initialization  from  noisy  data.  Resulting  model  for  (b)  isotropic  reconstruc¬ 
tion,  and  (c)  anisotropic  reconstruction  with  a  =  10. 


Figure  5.3:  Anisotropic  reconstruction. 


Figure  5.4:  Results  for  the  surface  area  reconstruction  with  weights  (a)  1,  and  (b)  10, 
isotropic  reconstruction  with  weights  (c)  1,  and  (d)  10,  anisotropic  reconstruction  with 
weights  (e)  1,  and  (f)  10. 


The  current  shortcoming  of  this  method  is  the  computational  speed  which  was  approxi¬ 
mately  one  hour  on  a  Intel  Xeon  1.7Ghz  Proc.  for  the  examples  presented. 


(c) 

Figure  5.5:  (a)  One  of  12  range  images  used  in  experiment,  (b)  result  of  feature  preserving 
reconstruction  from  a  similar  view  point,  and  (c)  a  different  view  point. 


Chapter  6 
Conclusion 


We  derive  a  variational  generalization  of  P&M  anisotropic  diffusion  for  feature  preserving 
surface  reconstruction.  This  generalization  is  based  on  a  robust  penalty  on  surface  normal 
vector  variations,  which  is  shown  to  have  important  advantages  over  using  surface  area 
and  the  quadratic  penalty  on  surface  normal  vector  variations  for  regularization.  The  data 
term  is  independent  of  the  prior,  the  ideas  introduced  in  this  paper  can  be  applied  to  other 
forms  of  surface  reconstruction  such  as  applications  in  tomography  [29].  We  use  implicit 
surfaces,  representing  the  implicit  function  on  a  discrete  grid,  modeling  the  deformation 
with  the  method  of  level  sets.  Therefore,  the  method  applies  equally  well  to  surfaces  that 
can  be  represented  in  a  volume.  The  results  shown  in  this  paper  are  not  possible  with 
previous  methods  in  the  literature. 

Measures  on  surface  normal  variations  require  solving  fourth-order  PDEs  on  level  sets. 
However,  by  processing  the  normals  separately  from  the  surface,  we  can  solve  a  pair  of 
second-order  equations  instead  of  a  fourth-order  equation.  This  method  is  numerically 
more  stable  and  computationally  less  expensive  than  solving  the  fourth-order  PDE  directly. 
The  shortcoming  of  this  method  is  the  computation  time;  however,  the  process  lends  itself 
to  parallelism,  and  therefore,  the  use  of  multi-threading. 
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